#include <iostream>
#include <cmath>
#include <arm_neon.h>
#include <semaphore.h>
#include <stdio.h>
#include <sys/time.h>
#include <algorithm>
#include <mpi.h>
#define N 4800
using namespace std;
int n=100;
float m[N][N],init_m[N][N];

void init(int n){
    for(int i=0;i<n;i++){
        for(int j=0;j<n;j++)
            init_m[i][j]=0;
        init_m[i][i]=1.0;
        for(int j=i+1;j<n;j++)
            init_m[i][j]=rand()%1000+1;
    }
    for(int k=0;k<n;k++)
        for(int i=k+1;i<n;i++)
            for(int j=0;j<n;j++)
                init_m[i][j]=int((init_m[i][j]+init_m[k][j]))%1000+1.0;
}
void reset(int n){
    for(int i=0;i<n;i++)
        for(int j=0;j<n;j++)
            m[i][j]=init_m[i][j];
}

void serial(int n){ //串行
    struct timeval h, t;
    float time = 0.0;
    gettimeofday(&h, NULL);
      for(int k=0;k<n;k++){
          for(int j=k+1;j<n;j++)
              m[k][j]/=m[k][k];
          m[k][k]=1.0;
          for(int i=k+1;i<n;i++){
              for(int j=k+1;j<n;j++)
                  m[i][j]-=m[i][k]*m[k][j];
              m[i][k]=0;
          }
      }
    gettimeofday(&t, NULL);
    time += 1000*(t.tv_sec - h.tv_sec) + 0.001*(t.tv_usec - h.tv_usec);
    cout<<"serial: "<<time<<" ms"<<endl;
    return;
}

void mpi()
{
    int thread_num;
    int rank;
    MPI_Comthread_num(MPI_COMM_WORLD,&thread_num);
    MPI_Comrank(MPI_COMM_WORLD,&rank);
    struct timeval h, t;
    float time = 0.0;
    MPI_Status status1;
    gettimeofday(&h, NULL);
    for (int k = 0; k < n; k++) {
        vector<int> b;
        vector<int> e;
        //初始化b和e
        for(int i=0;i<thread_num;i++) {
            b.push_back(0);
            e.push_back(0);
            int width = (n - k - 1) / thread_num;
            b[i] = k + width * i + 1;
            e[i] = k + (i + 1) * width + 1;
        }
        e[thread_num-1]+=(n - k - 1) % thread_num;
        for(int i=0;i<thread_num;i++){
            if ((n - k - 1) / thread_num == 0) {
                b[i] = (i == 0) ? k + 1 : 0;
                e[i] = (i == 0) ? n : 0;
            }
        }
        //除法并将结果发送给各个线程
        if(rank==0) {
            for (int j = k + 1; j < n; j++)
                m[k][j] /= m[k][k];
            m[k][k] = 1.0;

            for(int i=1;i<thread_num;i++){
                MPI_Send(&m[k][0],n,MPI_FLOAT,i,77,MPI_COMM_WORLD);
            }
        } else{
            MPI_Recv(&m[k][0],n,MPI_FLOAT,0,77,MPI_COMM_WORLD,&status1);
        }

        for (int i = b[rank]; i < e[rank]; i++) {
            for (int j = k + 1; j < n; j++)
                m[i][j] = m[i][j] - m[i][k] * m[k][j];
            m[i][k] = 0;
        }


        MPI_Status status2;
        int size=(e[rank]-b[rank])*n;
        if(rank!=0){
            MPI_Send(&m[b[rank]][0],size,MPI_FLOAT,0,99,MPI_COMM_WORLD);
        }else{
            for(int i=1;i<thread_num;i++){
                int i_size=(e[i]-b[i])*n;
                MPI_Recv(&m[b[i]][0],i_size,MPI_FLOAT,i,99,MPI_COMM_WORLD,&status2);
            }
        }
        MPI_Barrier(MPI_COMM_WORLD);

        MPI_Bcast(&m[0][0],n*n,MPI_FLOAT,0,MPI_COMM_WORLD);
    }
    gettimeofday(&t, NULL);
    time += 1000*(t.tv_sec - h.tv_sec) + 0.001*(t.tv_usec - h.tv_usec);
    cout<<"mpi: "<<time<<" ms"<<endl;
}


int main(){
    cout<<"n="<<n<<endl;
    init(N);
    reset(n);
    serial(n);
    reset(n);
    mpi();
    n=200;
    cout<<"n="<<n<<endl;
    reset(n);
    serial(n);
    reset(n);
    mpi();
    n=400;
    cout<<"n="<<n<<endl;
    reset(n);
    serial(n);
    reset(n);
    mpi();
    n=800;
    cout<<"n="<<n<<endl;
    reset(n);
    serial(n);
    reset(n);
    mpi();
    n=1600;
    cout<<"n="<<n<<endl;
    reset(n);
    serial(n);
    reset(n);
    mpi();
    n=3200;
    cout<<"n="<<n<<endl;
    reset(n);
    serial(n);
    reset(n);
    mpi();
    n=4800;
    cout<<"n="<<n<<endl;
    reset(n);
    serial(n);
    reset(n);
    mpi();
    return 0;
}
